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In this paper we derive the probability of the radial profiles of spherically symmetric inhomo- 
geneities in order to provide an improved estimation of the number density of primordial black holes 
(PBHs). We demonstrate that the probability of PBH formation depends sensitively on the radial 
profile of the initial configuration. We do this by characterising this profile with two parameters 
chosen heuristically: the amplitude of the inhomogeneity and the second radial derivative, both 
evaluated at the centre of the configuration. We calculate the joint probability of initial cosmologi- 
cal inhomogeneities as a function of these two parameters and then find a correspondence between 
these parameters and those used in numerical computations of PBH formation. Finally, we extend 
our heuristic study to evaluate the probability of PBH formation taking into account for the first 
time the radial profile of curvature inhomogeneities. 



I. INTRODUCTION 

The idea that large amplitude matter overdensities in the universe could have collapsed through self- 
gravity to form primordial black holes (PBHs) was first put forward by Zeldovich and Novikov l], and 
then independently by Hawking [2], more than three decades ago. This theory suggests that large amplitude 
inhomogeneities in the very early universe overcome internal pressure forces and collapse to form black holes. 
A lower threshold for the amplitude of such inhomogeneities <5 t h = (<Wp)th, was first provided by Carr 0], 
giving S t h ~ 1/3 at the time of radiation domination. This value was found by comparing the Jeans length 
of the overdensity with the scale of the cosmological horizon at the time of formation. 

The mass fraction of the universe turning into PBHs of mass M at their formation time, /?pbh(-^0, is 
computed using the probability density function (PDF) of the relevant field of inhomogeneities, which is 
provided by the cosmological theory The mass fraction /3pbh(-^0 is customarily given by the integral of 
this PDF over the amplitude 5 = Sp/p, with a lower limit equal to 5th [1, 0]- 

The probability of PBH formation is a useful tool to constrain the mean amplitude of inhomogeneities 
on scales which cannot be probed by any other methods. The PBH contribution to the energy density 
increases with time during the radiation-dominated epoch. For this reason, the PBHs formed considerably 
before the end of radiation domination are the most relevant to cosmology [B|, H, 0, H, H, EH • We will focus 
our study on these kind of PBHs and assume that the background matter at the time of PBH formation is 
radiation-dominated. To make this cosmological tool more precise, we must improve the calculations of the 
probability of PBH formation. This demands a more accurate evaluation of the threshold value 5 t h, or the 
equivalent curvature inhomogeneity 7?.th [HI, E2, EH, EH, El, EH ■ In search of these values, it was evident that 
the process of PBH formation depends on the pressure gradients in the collapsing configuration in addition 
to the amplitude [H, E3, EH ■ It was also found that such pressure gradients can modify the value of Sth 
significantly. Hence, when calculating the probability of PBH formation, one should consider the shape 
and radial profile of the initial configurations. These profiles are directly related with the internal pressure 
gradients. This is the main motivation for the present work. Previous studies, concerning the probability of 
PBH formation, take the amplitude of perturbations to be the only parameter determining the probability 
density. In addition to that, we include for the first time a parameter related to the slope of curvature profile 
at the edge of the configuration 1 . 



* ic.hidalgo@qmul.ac.uk Ja.g.polnarev@qmul.ac.uk 

1 In the context of dark matter haloes the question of initial profiles is effectively irrelevant because galaxies are formed from 
prcssureless configurations. The density profiles and shapes of virialized haloes result from the evolution of the initial high 
peaks and are not linked to the profile of initial configurations that we investigate here (see e.g. Ilij for a review on the 
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In this paper, we calculate the probability of finding a curvature configuration with a given radial profile. 
As follows from [U HI, Hi], [25 1 , PBH formation takes place only from nearly spherical configurations, so 
in this paper we restrict ourselves to the spherically symmetric case. In this first approximation we describe 
the radial profiles by introducing two parameters: the central amplitude of the curvature inhomogeneity 
lZ(r = 0) and the second radial derivative at the centre lZ"(r = 0), which is chosen mainly to avoid technical 
difficulties. The introduction of these parameters is a first step towards the full parametrisation of profiles 
in terms of even derivatives at the centre of configurations, i.e. in terms of 'RS 2n \0) (the odd derivatives 
7£(2n+i)^Q-j are a jj zer0 due to the assumed spherical symmetry). In the future, with the results from more 
accurate codes simulating the formation of PBHs, we will have at hand a larger number of conditions for the 
collapse of a curvature profile. An equal number of parameters will be required for the complete description 
of these profiles and the probability of finding them. In the meantime, only families of curvature profiles 
described by two parameters are available. We consequently limit ourselves to the two-parametric description 
of initial curvature profiles. 

The central amplitude 1Z(0) has been used in previous calculations of gravitational collapse 16], and the 
probability of PBH formation [H, [26|. In the present paper we compute the probability to find a given 
configuration as a function of the two parameters [1Z(0), 1Z"(0)]. We subsequently illustrate how this two- 
parametric probability is used to correct the probability of PBH formation. Such an exercise is presented 
for illustration purposes. The results, based on a non-rigorous but physically meaningful determination of 
the parameters which describe the initial profiles, show how the corrections to /?pbh are significant and they 
will be considered in more detail in future studies of PBH formation. 

This paper is organised as follows. In Section[TT]we calculate the joint probability distribution for 11(0) and 
1Z"(0). In Section UTIl we relate these parameters to those used in the most recent numerical computations of 
PBH formation. In Section lTvl we present the total probability of PBH formation, integrating the probability 
distribution derived in Section HT1 over the relevant region of parameter space [1Z(0) , 1Z"(0)]. In Section fVl 
we summarise our results and discuss future research in this area. 



II. PROBABILITY OF RADIAL PROFILE PARAMETERS OF COSMOLOGICAL 

PERTURBATIONS 

The most striking prediction of the theory of cosmological inflation is that initial quantum fluctuations 
are transfered into the inhomogeneities and structures observed in the universe today. After inflation, the 
universe is mostly flat with inhomogeneities of small amplitude on average (for a review, see [27]). Some of 
the inhomogeneous regions, present high amplitude (non-linear inhomogeneities) and these are the object of 
study in the present paper. Formally, the high amplitude inhomogeneous profiles describing configurations 
which collapse into PBHs are not perturbations. However, such regions are included in the statistics of 
random primordial curvature perturbations. That is, the statistics of random fields can be used to estimate 
the probability of finding high amplitude inhomogeneities. 

To describe large-amplitude inhomogeneities, we consider the non-linear curvature field lZ(t,r), as first 
described in [28[ . The non- linear curvature TZ(r, t) represents the relative expansion of a given local patch 
of the universe with respect to its neighbouring patches [2i|. It is described by the metric 

ds 2 = -N 2 (t, r) dt 2 + a 2 (t)e 2TC(t ' r) 7y (dr l + N*(t, r) dt)(dr J + N j (t, r) dt), (1) 

where a(t) is the scale factor and the gauge dependent functions N and N l are the lapse function and shift 
vector, respectively. These variables are determined by algebraic constraint equations in terms of the matter 
density p and pressure p and the metric variables 1Z, a and 7y. 

In this work we consider the non-linear configurations which correspond to non-zero and large 1Z inside 
some restricted volume, and zero outside, where the expansion of the universe follows the background 
Friedmann- Robertson- Walker (FRW) solution. There are several advantages of working with metric ((T|). 
First, 1Z is defined as a gauge-invariant combination of metric and matter variables (30| . Second, with the aid 
of the gradient expansion of the metric quantities [28, 31, 32, 33], lZ(r, t) is presented in the Einstein equations 
in a non-perturbative way. Finally, 1Z does not depend on time for scales larger than the cosmological horizon, 
as proved in [2^ . l3~ij and in a more general case in (35| . In the present paper we work in the superhorizon 
regime, where the field lZ(r) can be assumed to be time-independent. 



profiles of dark matter haloes and 20] for some recent results on this topic). 
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The primordial field of random perturbations we use follows a Gaussian probability distribution. In other 
words, we assume that the probability of finding a curvature perturbation 72 of mean amplitude i? is given 
by the probability density function (PDF) 



KfB *{-^)> (2) 
where T,-ji is the dispersion of the perturbation field 72, defined in terms of the two-point correlation function 

by 

Z 2 n(r H ) = (n(r)K(r)) = J A\nk W 2 (k, k H )V(k), (3) 

where W(k, ku) is the window function which smooths the field over spherical regions of size r# = 2ir /kn, 
the Hubble radius. The power spectrum of 72., V(k), is an output of the underlying cosmological model, as 
reviewed in [27| |. 

In more general cases, PDFs include the contribution of higher-order correlations (i.e. (727272) and all 
other correlations). To calculate such PDFs, a new formalism is required, such as that developed in [361 ]. 
Several studies have shown that the non-Gaussian correlations can sensibly mo dify the PDF of the amplitude 
of perturbations and consequently the number count of astrophysical objects [13, [H, [39| and PBHs fiol. l4ll . 
52, S3| when large non-Gaussianities arise in the primordial field of curvature fluctuations [HI, E^, I46L l47l]. 
Here however, we restrict ourselves to the Gaussian case where the PDF presents the form of Eq. @. 

In the following we calculate the joint probability of finding an amplitude 72(0) and the second derivative 



K"(0) 



.2 



(4) 

r=0 



at the centre, using the method developed in [36J . In order to compute the probability of a specific property 
of 72(r), we integrate the original PDF, which encodes all the information about the field, with the Dirac 
(5-functions of relevant arguments. In particular, the probability that 72.(0) = i?o is given by, 

P(tfo) = / [dn]¥{n) 6 [72(0) - tf ] , (5) 

were [dlZ] indicates integration over all possible configurations 72(k) in Fourier space. Hereafter we consider 
72(0) and 72." (0) as statistically independent parameters. Hence, the probability of having 72." (0) = #2, is 
given by the integral 

,)= [[ctll]F(1l)6[1l"(0)-# 2 ], (6) 

In the rest of this section we show roughly how this method works. The details of the following results are 
presented in appendices O and [B] First we expand the smoothed curvature perturbation profile 72.(r) in 
terms of spherical harmonic functions: 

/H 3 k 
— 72(k)exp(ik.r), (7) 

with 

00 i 00 

^( k )=E E J2 n "n Y im(0,<P)Mk)- (8) 
i=0 m=-i n=l 

Here Yg m are the usual spherical harmonics on the unit 2-sphere and ip n (k) are a complete and orthogonal 
set of functions in an arbitrary finite interval < k < A (for an explicit expression of ip(k) see appendix 
IX]) . It is worth mentioning that the cutoff A is imposed to artificially compactify the momentum space. This 
allows us to provide an explicit definition of the functions ipn(k) and a complete set of functions ip for the 
expansion of 72 (k). In turn this condition allows a regularisation of the path integral J[dR] by considering 
the harmonic expansion © in a finite interval in Fourier-space < k < A. At the end of the calculation 
we can take the limit A — > 00 and the results will remain unchanged. The coefficients in the expansion are 



generically complex, so we separate real and imaginary part introducing TV^ n = a™ n + ib'^ n . The reality 
condition for the curvature field, 7?.(k) = 1Z(— k), is met when 



— m / 1 \^+m m 



(-1) 



^|n- 



(9) 
(10) 



In particular, the m = modes require a°| n and &°| n to be zero for odd and even £, respectively. After 
evaluating the expansion (0-© at 7Z(r = 0), we can use the relation, 



dVLYJ 



where dfl = sm(8)ddd<p, and integrate Eq. to obtain 



47T^ 



mOi 



(11) 



ft(0) 



3 £- /" k 2 dk °° £ °° / 

exp (ik . r) , r=0 = / ^ E E E ^ (V^o 



(2 



>m0 



=E«o 

71=1 



(12) 



To evaluate the central second derivative, we follow the same steps and obtain 

d 3 k 



(27T)' 



r ^(fc)(ifc) 2 exp(ik-r)| r= o 



dk 



V^(2tt) 



r ^„(fc)fc 4 = ^ 2 



(13) 



The intermediate steps of the derivation of Eq. (|13[) are presented in Appendix [XJ 

To proceed with the computation of the probabilities given in Eqs. ([5]) and ^ we must integrate over 
all configurations in Fourier space. With the aid of the expansion © we can express the measure of such 
integral in terms of the expansion coefficients satisfying the reality conditions and (JTTJJ) , i.e., for any ^[72.], 
functional of lZ(k), the following integral can be represented as 



V[n] [dTZ] 



oo £ oc 

nnn^ 

£=0 m=l n=l 



9[R] da} 



OO OC pQQ /»CC 

n n * / d °u / d6 2 P +ii 9 

„_n „_ i «^ — OC J — OO 



p=0 g=l 



(14) 



where the constants /i and /i are weight factors. In our calculation of probabilities, such factors are absorbed 
by the final normalisation of the joint probability. 

As mentioned before, we restrict ourselves to the Gaussian PDF. In terms of the spherical harmonic 
coefficients (see Appendix |B|) . this means that 



( 1 OO i OO OO OO 

E E E + i^J 2 - E E la °^ 
v ' i=0m=0n=l y ' p=0 q=l 



\ 2 + K +1 J 



(15) 



In order to obtain the probabilities of the mentioned parameters from Eqs. ([5]) and ©, we use the repre- 
sentation of the Dirac (5-function 



6(x) 



/OO 
dz exp[iz x] 
-co 



(16) 



This allows us to write, for example, the (5-function in Eq. (|5|) in terms of the spherical harmonic coefficients 
as 



S(TI(0) -0 O ) = / dzexp 



E4 



dkip n k — 



2 (2^ 3 



4tt 



(17) 
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In the same way, the representation of 6 [72." (0) — can be written in terms of harmonic coefficients with 
the aid of Eq. IQ3J). 

We now have all the elements needed to derive the probability of the parameters 11(0) and 72"(0). Substi- 
tuting expressions (Tl5|) and (fT7|) in Eq. ((5|) , we perform the functional integral with the aid of the decompo- 
sition (| 14[) . In this process we discard all the Gaussian integrals because they contribute to the probability 
only with a multiplicative constant which will be included in the final normalisation. On the other hand, the 
Dirac 5-function contributes with exponential factors of a^ n to the integrals. The integrals of these param- 
eters are computed by completing squares of the exponential arguments, so the integrals of such coefficients 
include a set of shifted Gaussian functions (see Appendix [B] for the details of this procedure) . The integral 
(lU) can be performed following the same steps and using the corresponding expressions (fT3]) , (fT4|) and (fT5|) . 
The final probability density for the pair of parameters 11(0) and 72." (0), is the product of the integrals (fT2|) 
and (TT3j) . i.e. 

P (72(0) = d , 11" (0) = $ 2 ) = A exp ^-^f - - , (18) 

where £( 2 ) an d £(4) are the dispersion of the amplitude and the second derivative respectively, and A is a 
normalisation factor obtained from the condition that the integral of the joint PDF over all possible values 
of the two independent parameters equals unity. The final normalised joint probability density is 

**•*> = ^m*m^ ("4; (19) 

According to the Press-Schechter formalism of structure formation [4( , the PDF is integrated over all per- 
turbations which collapse to form the astrophysical objects under consideration. In this way we calculate 
the mass fraction of the universe in the form of such objects. To apply this formalism and calculate the 
probability of PBH formation and integrate the PDF (|19[) , we require the range of values 72(0) and 72" (0) 
which correspond to PBH formation. In the next section we will obtain this range with the help of the 
results of numerical computations presented in (l8l | . 



III. THE LINK BETWEEN PERTURBATION PARAMETERS AND THE CURVATURE 
PROFILES USED IN NUMERICAL CALCULATIONS 



A. Initial conditions 

As demonstrated by the first numerical simulations of PBH formation [l7]], whether or not an initial 
configuration with given curvature profile leads to PBH formation, predominantly depends on the following 
two factors: 

1) The ratio of the size of the initial configuration ro to the size of the closed universe — 

a(t) L dr/yl — r 2 (evaluated at the initial time), which is a measure of the strength of gravitational field 
within the configuration. 

2) The smoothness of the transition from the region of high curvature to the spatially flat FRW universe, 
which is characterised by the width of the transition region at the edge of the initial configuration and it is 
inversely proportional to the pressure gradients there. Strong pressure gradients inhibit PBH formation. 

The numerical computations presented in [l8| (hereafter PM) give the time evolution of the configurations 
with initial curvature profiles accounting for the above-mentioned factors and collapsing in a radiation- 
dominated universe. In that paper the initial conditions are obtained with the help of the quasi-homogeneous 
asymptotic solution valid in the limit t — > 0. This solution to the Einstein equations was first introduced 
by Lifshitz and Khalatnikov [35[ (see also [H, HH). Following [T^], PM used this asymptotic solution to 
set self-consistent initial conditions for curvature inhomogeneities, the initial curvature inhomogeneity being 
described by the spherically symmetric curvature profile K(f). This sets the initial conditions for the process 
of black hole formation. Asymptotically, the metric can be presented in terms of K(f) as 

ds 2 = s 2 (n) l-drf + ^-^df 2 + f 2 \d9 2 + sin 2 6d(j) 2 ] I , (20) 



G 



where s(rj) is the scale factor, rj is the conformal time and we write f for the radial coordinate to distinguish 
it from the coordinate of the metric |T]) . An advantage of working with this metric is that it contains the 
curvature profile K(r) explicitly. We choose a set of coordinates with the origin at the centre of spherical 
symmetry and fix K(0) = 1. The condition that K(f) is a local inhomogeneity requires that K(f) = for 
radii f larger than the scale fo, where the metric matches the homogeneous FRW background. 

In PM the profiles K(f) are presented in two forms, one of which is characterised by two independent 
parameters a and A as 



K(f) = 



1 



2A 2 



exp 



2A 2 



(21) 



The results of the numerical simulations in PM indicate that PBHs are formed in the region of the parameter 
space [a, A] shown in Fig. [T]a,. 
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FIG. 1: a) The left plot shows the parameter values for initial configurations which collapse to form black holes. A 
characterises the width of the Gaussian curvature profile, while a characterises the deviations from a Gaussian profile, 
as can be seen in Eq. JS). b) In the [K(Q), TZ"(0)] plane three regions of integration are considered to compute the 
probability of PBH formation. Region I is the region enclosed by the solid curves and corresponds to the region noted 
by BH in Fig. [lji. Region II is the region to the right of the grey dotted line re pres enting the surface of integration 
considered in previous studies where only the amplitude is taken into account [2g|. Region III is the region above 
the solid line and between the dashed lines. The physical characteristics of profiles with values in this region are 
described in section IIII CI 



B. Physical criteria for the identification of parameters 



We proceed by finding the correspondence between the two sets of parameters, [72.(0), TZ" (0)] and [a, A], 
both of which describe the initial curvature profiles. Assuming that the size of the configuration, tq, is 
much larger than the Hubble horizon r# = H , where H is the Hubble parameter, we can use the gradient 
expansion of the functions in metrics (fl]) and (|20p . In this case, the time derivative of any function f(t,r) 
is of order f /t ~ H / and significantly exceeds the spatial gradient which is of order / /tq. Hence the small 
parameter in the gradient expansion is 

«=? = 4> (22) 
?*o all 

where k is the wave-number corresponding to the scale of the configuration. Taking into account that e — > 
when t — > 0, one sees that the gradient expansion is very similar to the quasi-homogeneous solution (35j. 
For the metric (JTJ, using the coordinate freedom to set N l = and ignoring any tensor contributions, i.e., 



taking 7^ = <5y, the expansion of the Einstein equation G ° — 8irGT Q ° to order e 2 can be written as 2 , 



I 1^ _{s) R _^_ {N _ 1) \ +0(e 4 )=8KG(po + 6p)+0(e% 



2 V a 2 



(23) 



where ^ 3 'R is the spatial curvature, or the Ricci scalar for the spatial metric gij. To order zero in e, we have 

3d 2 



8irGp , 



(24) 



which corresponds to the homogeneous part of (|23|) . As shown in [lfj, l5l|, [52j , the time slicing can be set to 
a uniform expansion gauge in which 



(25) 



where T — 1 is the sound-speed squared. Using (|23|) . (f24|) and (|25l) . we find the equivalence between the 
spatial curvature and the matter overdensity 



3 H V 3r 



In consequence, the gradients establish a correspondence with the pressure gradients 



V(3)R= 8^4 + 3T V(M 8*0 / 4 + 3F 



3r 



3 V 3 r(r-i) 



(26) 



(27) 



where V = {g rr )~ x / 2 d/dr. The last equation shows that the gradient of the spatial curvature is directly 
related to the pressure gradient. Hence, subject to these two physical conditions at the edge of the config- 
uration, we relate the profiles TZ(r) and K(f) by equating the spatial curvature and its gradient for metrics 
© and d20D- That is, 



2U"{r) + {K'{r)f exp(-2ft(r)) = 3K(f) + rK'(r), 



(28) 



and 



I d 



(( 3 )R) 



[B'R" + ft"'] exp(-3ft(r)) = 



"I - Kr 2 " 


1/2 ^ 







2fK'(f) + -r 2 K"(f) 



(29) 



By definition of the edge of curvature configuration, the three curvature must vanish at this point, so Eq. 
implies 



2K"(r ) + (W(r Q )) 2 = 0, 



and 



3K(f )+f K'(fo) = 0. 
As a consequence of this, the gradient relation (j29|) can be written as 

1/2 



[W(r Q ) 3 - 271"' (ro)] exp(-3ft(r )) 



I-Kf 2 



-12K(f )+f 2 K"(fo)]. 



(30) 
(31) 

(32) 



2 For the complete second order expansion of the metric quantities, see for example, |29l l32| . 



This establishes a relation between 72(r) and K(f) at the edge points r$ and r$. The configuration K(r) is 
parameterised by [a, A], as shown in Eq. (|2ip . As follows from condition ([31]) (see also PM), the radius ro 
can be written in terms of those parameters as 



5a- 2 + ^J{ba - 2) 2 - 24a 



2a 



A 2 . (33) 



Then we use two more equations obtained from the conformal transformation of coordinates at zero order 
in e: 

aW TCM ^ S 2 (,) T -f^ (34) 

and 

a 2 (r) e 2TC « r 2 dtf = S 2 ( V ) f 2 dQ 2 . (35) 

Because the homogeneous Einstein equations are identical in both metrics, the scale factors a(r) and s(r;) 
can be identified, a(r) = s(rj). Thus we find a relation between the radial coordinates, 

e n{r) r = f (36) 

and an integral relation between the configurations 



e 



(i 



n{x) dx = r dx (37) 
Jo y/l-K(x)x 2 ' 



One can verify that Eqs. (|28p. (|2^)) and (|3"T|) are not independent. For example, Eq. (|2U)) follows from 
and f37]) . 

In the previous section we have developed a method to account for the probability of any set of parameters 
describing the curvature profile. For simplicity we have chosen the pair [72.(0), 72." (0)]. We now illustrate 
how to relate [72(0), 72" (0)] and [a, A] by considering the parabolic profiles 

72(r) =72(0) + i72"(0)r 2 . (38) 

This parametrisation meets the minimal requirement of covering the [a, A] parameter space in Fig.[T£i. 
Eqs. ((30|) . (I37|) and (|36|) are now reduced to the following system of algebraic equations: 

r ° = ^T^O)' (39) 
m - 21og [^IWJ-V (1 _ K ^ 2)1/2 ) , (40) 

R = (41) 

where fp is given in terms of [a, A] by Eq. (|33[) . 



C. Parameter values leading to PBH formation 

As follows from the numerical computations [IH which used the parametrisation (|2"Tj) , PBHs are formed 
in the [a, A] region shown in Fig. QJu Equations (|40|) and (|4Tj) map this region to the Region I in the 
space of parameters [72(0), 72." (0)] shown in Fig. [TJd. The Jacobian of the transformation corresponding to 
this mapping is non-vanishing, which guarantees a one-to-one correspondence of the region 'BHs' plotted 
in Fig. with Region I in Fig. []J>. Each point here corresponds to a parabolic profile which leads to the 
formation of PBH. 

For each one of these parabolic profiles, there is a family of non-parabolic profiles with the same central 
amplitude 72(0), the same configuration size r , and the same behaviour near the edge, as shown in Fig.O 
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FIG. 2: The curvature profile for three different families of configurations with common central amplitude 1Z(0) — 1. 
The configurations shown by the dashed lines have value of 1Z" (0) larger in absolute magnitude than the parabolic 
one shown in black. The configurations shown by the dotted lines have a value of TZ" (0) smaller than the parabolic 
one. All profiles satisfy conditions (|30[l and (I32|l . 
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FIG. 3: The logarithmic probability of PBHs (3 calculated using Eq. (42) with a power spectrum with two tilt values 
(n s = 1.32 for the left plot, n a = 1.47 for the right plot). The lines show the integration for the three different regions 
sketched in Fig. f. The integral over the region I (/3i(M)) corresponds to the dashed lines, and the integration over 
the region II (/3n(M)) to the solid lines. The probability integrated over the region III (/3m (M)) is represented by 
the dotted lines in both figures. 



In that figure, the profiles lying below the parabola correspond to larger absolute magnitudes of 1Z"(0) and 
do not form PBHs because they have lower average gravitational field strength and higher average pressure 
gradient. The non-parabolic profiles which lie above the parabolic one (with smaller absolute magnitude 
1Z"(0)) should also collapse to form PBHs because they correspond to higher average gravitational field 
strength and lower pressure gradient. 

In the parameter-space [71(0), 7Z"(0)], this last set of profiles corresponds to Region III in Fig. [TJd. This 
region will be included in the calculation of the probability of PBH formation in the next section. 



IV. TWO PARAMETRIC PROBABILITY OF PBH FORMATION 
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To calculate the probability of PBH formation, which is equivalent to the mass fraction of the universe 
going to PBHs of given mass, it is customary to use the standard Press-Schechter formalism J4[. This has been 
widely used in previous calculations of the one parametric probability of PBH formation pi IH. Icl l53l l54j . 
When the probability depends on a single amplitude parameter, this method reduces to the integration of 
the corresponding PDF over the relevant perturbation amplitudes. The final integral is equivalent to the 
mass fraction of PBHs of mass M ~ (T - lf/ 2 M H rs (T - l) 3 ^k M /(2n) 0, with the soundspeed x/T - 1 
measured at the time formation 3 . Here we extend the standard Press-Schechter formalism to derive a two 
parametric probability introducing the second derivative at the centre of the configuration as an additional 
parameter. When the [72/(0), K"(0)] region is a square [72i < 11(0) < K 2 , K'{ < 72." (0) < 72' 2 '], the integrated 
two parametric probability for objects of mass M is 

I3pbh(M) = ["* dd I" 12 <M 2 P(tf ,tf2) = 
J Tlx Jn'{ 



erf ( — = I — erf 



V2S (2) (Af)/ \V2S( 2 )(M) 



erf f *S -erf' ** 



(42) 



We use this result to integrate numerically over a mesh of small squares covering each one of the regions 
of the plane [72.(0), 72," (0)] shown in Fig. [1}d. We call the integral over region I (3i(M), and correspondingly 
the integrals over regions II and III are called (3\i(M) and f}\\\(M). The mass dependence of these betas 
for two different power-law spectra Vjz(k) cx fc" -1 , are shown in Fig. [3] As dictated by the Press-Shechter 
formulation, such integration corresponds to the fraction of mass density in the universe that has collapsed 
into objects with mass M. We remind the reader that our choice of n s and the mean amplitude is for pure 
illustration purposes. With the values used in this paper, copious amounts of black holes are produced. A 
red spectral index and power spectrum inferred from the CMB data corresponds to a low number of PBHs. 
However, to assume that the same values of the power spectrum and spectral index are valid on scales 
relevant to PBH formation (30 decades of mass below the mass scales correspondent to CMB observations) 
is a very strong extrapolation. At the present time we cannot exclude that the values of n s and the power 
spectrum are different than those given by the CMB. To explore the structure formation models that match 
CMB observation values and also produce considerable number of PBHs is a great task beyond the scope of 
our paper. This important issue is currently under investigation [lOL l5a | . 

We contrast the case of parabolic profiles described by Eq. ([55]) with the non-parabolic set presented in 
Fig. H]by plotting the ratios of probabilities /3i//?n and @m] fli for different values of V-r,- This is presented 
in Fig. [4j This figure shows that the probability of PBH formation can be larger than the previous one- 
parameter probability computed from the integration of Region II as done in previous studies (26|. This 
important result requires confirmation from more detailed numerical simulations of PBH formation in this 
region of parameter-space. The uncertainty is explained by the fact that the two parametric calculation of 
the probability of PBH formation is still incomplete. This should be complemented in the future by the 
introduction of all relevant higher derivative parameters and the higher-order correlations in the PDF. 



V. DISCUSSION 



We have developed a method for calculating the two-parametric probability of PBH formation, taking into 
account the radial profiles of non-linear curvature cosmological inhomogeneities. This is the fist step towards 
calculating the TV-parametric probability which takes into account the radial profiles more precisely than 
studies using the amplitude as the only relevant parameter. Using the results of sophisticated numerical 
computations, we obtain the values of 72" (0) that are relevant for PBH formation. Subsequently we have 
incorporated these values to the total probability of PBH formation. Finally we have provided an example 
of the consequences of this probability for the statistics of PBHs. 



Throughout this paper we consider configurations that collapse in a uniform radiation dominated background. Thus we use 
the value T = 4/3 



11 



5 1 ■ ■ ■ ■ r 




-2.5 -2 -1.5 -1 -0.5 



log 10 (Pn) 

FIG. 4: The horizontal axis of the figure is the amplitude of the power spectrum at scales relevant to PBH formation. 
The grey dashed line shows the ratio /3i//3ii where /3i is the probability density integrated over region I in the 
[TZ(0), 7?,"(0)]-parameter space of Fig. lb, and /3n is the corresponding probability integrated over region II of the 
same figure. The black line is the ratio /3m over /3u. The dotted line shows the reference case of the single-parameter 
probability. 



The results obtained show that, if we restrict ourselves to the PBH formation calculated for parabolic 
profiles only (as described in Section iLTlj) . then the total probability of PBH formation is orders of magnitude 
below previous estimations. On the other hand, with the aid of heuristic arguments we show that a much 
larger region of parameter-space [7^(0), TZ" (0)] representing non-parabolic profiles should also be considered 
in the estimation of the probability of PBH formation (see Fig. [2J. In this case, the total probability 
of PBH formation is higher than the single-parameter estimate of previous works. In this case, we have 
an opportunity to impose new bounds on the power spectrum on the scales relevant for PBH formation. 
Analysing the uncertainty of our results, mostly due to the heuristic nature of the present study, we have 
demonstrated how much we still have to understand about the formation and statistics of PBHs. The 
physical arguments supporting our results should be made rigorous by direct verification with numerical 
hydrodynamical simulations of PBH formation. This in turn would provide valuable support for the initial 
motivation of this work. 

The main conclusion of this paper is that the amplitude of initial inhomogeneities is not the only parameter 
which determines the probability of PBH formation. The ultimate solution of the problem requires a greater 
set of parameters and a larger range of their values to determine all high curvature configurations that form 
PBHs, which is a huge task for future research. In the meantime, we have a method to operate with the 
statistics of all these parameters. 
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APPENDIX A: HARMONIC DECOMPOSITION OF TZ AND FOURIER REPRESENTATION OF 

11(0) AND K"(Q) 



The Fourier expansion of the smoothed curvature perturbation profile TZ in terms of spherical harmonic 
functions is, 

oo i oo 

K(Kt)=J2 E Y. n T\nit)yUe k Ak)Mk)- (Ai) 

1=0 m =-£ n=l 

The radial functions of the harmonic decomposition can be defined by 



J„+i(a«) Afc 2 V "A 

where W(fc; /cm) is the window function with smoothing scale /cm and a™ is the n-th root of the Bessel 
function of order v, J v (k). Note that the functions of the radial coordinate in the expansion are the Bessel 
functions up to a factor. The set of functions ip n (k) is orthonormal under the product 

V(k)W 2 (k; k M ) 

and the completeness relation can be written with the aid of the Kernel in the internal product: 



h () 



V(k )W 2 (k M ) 



J2^n( k o)Mk) = S(k-k ). (A4) 



In Eq. (|A3|) A represents an artificial compactification of the momentum space which is used only to have 
an explicit definition of ip at hand. With the definitions above, the expansion (|Alj) describes the curvature 
perturbation with power spectrum V(k) smoothed over a scale /cm- 

Let us now construct a parameter to represent the second radial derivative of the field TZ in Fourier space. 
The first radial derivative of TZ(r) is 

JU(r) = Jd 3 kY, Kn^riOk, fe )Vn(fc)e i(k - r) X |-(ik • r). (A5) 

If we write the spherical coordinates in a Cartesian basis, we have 

k r = \k\\r\ {sm(9 r ) cos(</V) sin(#fc) cos(cj)k) + sin(0 r ) sm(cj) r ) sin(9k) sin(0fc) + cos(9 r ) cos(0k)} , (A6) 

with (4> r ,9 r ) the set of angles of the vector r and Ok) the corresponding pair for k. 

The scalar product (|A6[) can be simplified if we note that the integral in Fourier space spans all possible 
directions of k, so we can choose an arbitrary direction for r. In particular, we can fix the r so that 
cos(# r ) = 1. This allows us to write a simple expression for the radial derivative, 

^(kT) = |fc|cos(0 fe ). (A7) 

The first derivative of the profile TZ(r) at the centre of the configuration is zero by the symmetry of spherical 
configurations. By construction of the spherically symmetric Fourier modes, this condition is satisfied identi- 
cally at r = 0. The first non- vanishing parameter that gives information about of the profile of perturbations 
is the second derivative. Fixing the direction of the vector r in the scalar product we have 

g-^Tlir) = J d 3 k TZZJr{e k Ak)Mk) x ^cos^K^. (A8) 
With the standard definition of the spherical harmonics, 
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where the normalisation factor is used for orthonormality purposes, one can write the factor cos 2 (9k) as the 
sum of two spherical harmonics, 



cos 2 (6> fe ) = i\/47r 



1 



47T 



(A10) 



where the Y indicates the complex conjugate. 

Using the normalisation rule for spherical harmonics, 



dnY l m (e,<f>)Y k n (e,<f > ) = s mn s u , 



we can integrate Eq. (|A10[) in the derivative (|A8|) and arrive at the expression 



n"(o) = J dkk 4 



4ir 



(All) 



(A12) 



This is the result used in Section ITT1 

To finish this appendix we show how the introduction of a parameter off the centre, say 1Z' (ro), generates 
a large set of constraints on the values of the coefficients TZ^ n - The integral in Fourier space representing 
this derivative is 



K'(r)\ r -_ 



d A k 
d 3 k 



\k\ lZ(k) cos(^fc) exp [i k r cos(6>fc)] 

(icos(6» fc )r fc) i 



\k\K(k)Y? 



E 

.s=0 



(A13) 



where we have expanded the exponential function in Taylor series. Each power of cos(#) can be expressed 
in terms of spherical harmonic functions YP, This means that the last integral consists of a series of terms 
of the form, 



dfl YJ l (6, <j>)Y?(6, <p)Yg(6, <j>). 



(A14) 



Each of these integrals is a Clebsch-Gordan coefficient. These steps are enough to show that, while a 
parameter IV (r = r ) might represent an improvement in the estimate of the final probability of PBH 
formation, it takes a long calculation to complete squares, add normalisation factors for each coefficient 
TZ% n , and arrive at a final expression like eq. ((T9|) . This goes beyond the goals of the present paper. 



APPENDIX B: THE PROBABILITY OF 11(0) AND TZ"(0) 



At any time i, the probability distribution Pt[R], is formally obtained through the inverse Fourier transform 
of Z t [rj] , a generating functional which can be expanded in terms of the n-point correlation functions [36| , 

Z t {r)]=expJ2- ■■ d 3 y 1 ---d 3 y n r 1 (y 1 )---r 1 (y n )(K(t,y 1 )---'Jl(t,y n )), (Bl) 

n=0 H ' 

Hence, up to an overall normalisation, 

F t [K] cx J [6.1]] exp ^-i J d 3 x K(x)r)(x)J Z t [rj\. (B2) 

If all correlation functions of three and more points are set to zero, then W[R,] cx G[7£]. Assuming this for 
the generating functional (|B 1[) . the expression to integrate in the Fourier space is, 



Ptfaft] =exp 
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The functional integral of this expression gives the probability of 72. To solve this integral we complete the 
square of 77 factors and make the finite field redefinition 



ry(k) ^ f/(k) = 77(k) + i(27r) 3 



R(k) 



(K(t,-k)K(t, -*.))" 



(B4) 



where the prime ' attached to (72(i, k)72(i, — k))' indicates that the momentum-conservation ^-function is 
omitted. The measure [dry] is invariant under this shift, giving J[drj\ = J[dr)], whereas P t [77; 7^.] can be split 
into an 72-dependent piece, which we call <G t [72], and a piece that depends only on fj but not 72, 



t [ m n] » G t [K] ex P (--y 1 2 fj ( k , 1 „ i k 2 ) in ( / . k .ma. k 2 1 > 



where Gt[72] is a Gaussian in 72, 



G f [72] =exp d 3 fc!d 3 fe 2 (72(t,k072(t,k 2 )) — 



^(k!)7e(k 2 ) 



(^.kij^-ki))' 



(B5) 



(B6) 



When we make the expansion of the fields 72. (k) in terms of the spherical harmonics as in Eq. (|A1[) and using 
the explicit expression for the two point correlation Eq. ([3]) we obtain, 



G[72] = exp 



r / dfU k-dk 



k 3 



1 



(27T) 3 27T 2 7>(fc)W 2 (/fc) 



X E E ^|n 1 ^i^x-x(^^,n, a («^)^»i(*)V'n a (fc)J, 

£l ,1711,711 ^2,W12,"2 / 



(B7) 



where the normalisation factor has been left aside and can be recovered by demanding the integral over all 
values to be equal to 1. 

The harmonics Yg m and ip n integrate out of this expression entirely, using the orthonormality relation (|A3j) 
and the spherical harmonic completeness relation (|A1 1|) . Moreover, after rewriting the a and b coefficients 
with m < in terms of the m > coefficients, we obtain 



00 00 



= «P - E E E W + Kf - E E W + l 6 °+l|»| 2 • CB8) 

\ V ' e=0m=ln=l * 1 t=0 n=l J 



which is the gaussian expression presented in (|T5j) . 

In order to find the probability for given values of the central amplitude 72(r = 0) = i)q, we integrate G[72] 
with the (^-function factor in (fT7|) . We introduce the Fourier representation of the (5-function to write, 



0) oc I [dTZ] ( dz G[72] exp 

J J —OO 



^2 a 0\n^r, 



(2nf 



4tt 



da 



(B9) 



where the functional measure is understood to be Eq. (|14[) . The final answer is obtained by integrating out 
z together with all of the a and b coefficients. In order to achieve this, it is necessary to decouple a®,, z and 

i9n from each other by successively completing the square in Oq| and z. Working with Oq| first, we find 



1 1 



exp 



: exp 



4tt 2 (2tt) 3 ^ 

K ' 71=1 



EKJ 2 +i*E a oi„^ 



1 1 



^3 E ( fl o|n " i2^ 2 (27r)^I]„) - (2^) Vz 2 £ 2 2) 



4tt 2 (2tt) 3 ^ 



(BIO) 



where we have introduced a function £ 2 2 ), defined by E 2 ^ = X^^Li ^n- ln the final probability distribution, 
E 2 ^ will turn out to be the variance of 72(0). From Eq. (|B10|) . it is clear that making the transformation 
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l®, + i27r 2 (27r) 3 zS„ suffices to decouple a^ n from z. The measure, Eq. (H3J), is formally invariant 
under this transformation. We can also complete squares for the variables z and $q, giving 



exp ( -(2^) 3 7rVE 2 2) -i^)l^ z 



4tt 



exp 



-(2tt)VS 2 z + i ^ 

V 7 (2) I 2vT 2 \/4^I] 2 



l) 2 

2E 2 



(Bll) 



As before, the finite shift zhz - i$o/27r 2 \/47rI] 2 2 - ) leaves the measure intact and decouples z and $o- The 
a, 6 and z integrals can be done independently, but since they do not involve "do they contribute only an 
irrelevant normalisation to P(i9q)- The result is the Gaussian distribution for $0j 



P(i?o) oc exp 



2S (2) , 



(B12) 



It remains to evaluate the variance E 2 ^. In the present case, we have E„ = J A dfc k 2 ip n (k). From the 
completeness relation Eq. (|A4[) . it follows that 



(B13) 



£ fc v„(fc )fcv„(fc) = fc2p( y (fco) ^fc - **). 



S 2 2 j is now obtained by integrating term-by-term under the summation. The result coincides with the 
conventional smoothed variance, 



Jo 



dlnfc W 2 (Jfe;fcff)P(fc). 



(B14) 



Thus, as expected, Eq. (|B12[) reproduces the Gaussian distribution © which was derived on the basis of the 
central limit theorem, with the proviso that the parameters (such as E 2 2 p describing the distribution of $o 
are associated with the smoothed field 1Z. 

For the case of the central second derivative we integrate this probability against the 5-function containing 
the desired condition He 



P(0 2 ) = / [dR]G[K]6 [K"(0) - tf 2 



(B15) 



Using again the expression of the ^-function as an integral and the condition on the derivative in terms of 
the spherical harmonic coefficients we have: 



P(t? 2 )oc j [dTZ] J dzG[R] exp 



4tt 



r a 2\n + °0|n 



where the factor E^ 4 ' is defined as 



In the integral (|B16|) the terms with factors of Oq. are 

1 x v 

V ' n=l n=l 



(B16) 



(B17) 



exp 



(4) 



(B18) 



Completing squares, this last expression is equal to 



exp 



-i^2CTE(l a o|J-i(2-) 3 2 7 r 2 zE( 

^ ' n— 1 



- (2^) We 



(4) 



(B19) 
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In the same way we can complete squares for the expansion factors a^ n 

n—l n—1 

1 00 / A 2 



exp 



'47r 2 (27r) s 



= exp 



(B20) 



And finally one can also complete squares for the terms containing the variable z which are independent of 

a 0\n and a 2|„' 

.3(2tt) 3 



exp 



exp -(2-7r) 3 7r 2 
where for simplification we have written 



47T 

5 tfa 



5 1? 2 



(4). 



2 y 2 



(B21) 



J (4) 



1(4) 



(B22) 



To evaluate this variance of the second derivative we use again the property (|B13[) to integrate the complete 
summation and obtain 



^(4) 



So by making the change of variables 



and 



d\nkW 2 {k,k H )V(k) k 4 



4tt 2 

a§, n ~a°|„ + i^(2,r) 3 £W* 
.5 d 2 



z 1— >z + 1 



12^ S 2 



(B23) 

(B24) 
(B25) 

(B26) 



(4) 



we can perform all the integrals and eliminate the gaussian ones which contribute only up to an overall 
numerical factor subsequently absorbed by normalisation. The remaining factor expresses the probability of 
finding a perturbation 1Z with a central second derivative of value #2 > 



'[Tl"(r = 0) = 2 ] oc exp 



Ml 
2S( 4) 



(B27) 
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